rm(list = ls(all = TRUE))
graphics.off()
library(MASS)
library(Matrix)
library(lubridate)
set.seed(07031969)
options(warn=2)

load("ONScars.RData")
sel <- which(data[,"year"]>=88)
data <- data[sel,]

#######################
# Price distributions #
#######################
price <- data[,"price"]/1000
from <- 1000/1000
to <- 20000/1000

year 		<- data[,"year"]
month 		<- data[,"month"]
quarter 	<- floor((month-1)/3)+1
yrq 		<- year*100 + quarter
quarters 	<- sort(unique(yrq))
T 			<- length(unique(yrq))

fx <- NULL
average.price <- NULL
quartz(width=9,height=7)
pdf(file = paste("..//Manuscript//Figures//Figure_1a.pdf",sep=""), onefile = TRUE,width=7,height=7)  										# FIGURE 1(a)
for (t in 1:T){
	sel <- which(yrq==quarters[t])
	fx <- cbind(fx,density(price[sel],from=from,to=to,n=100,bw=0.75)$y)
	xvals <- density(price[which(!is.na(price[sel]))],from=from,to=to,n=100,bw=0.5)$x
	colour <- (paste("grey",100-(40+2*t),sep=""))
	colour <- (paste("grey",(3*t),sep=""))
	if (t==1){
		plot(xvals,fx[,t],typ="l",xlab="£ '000's",ylab="Density",col=colour) 
	} else {
		lines(xvals,fx[,t],col=colour)
	}
	average.price <- c(average.price,mean(price[sel]))
}

dev.off()
graphics.off()

quartz(width=9,height=7)
pdf(file = paste("..//Manuscript//Figures//Figure_1b.pdf",sep=""), onefile = TRUE,width=7,height=7)											# FIGURE 1(b)
contour(seq(88.125,95.5,0.25),xvals,t(fx*10000),ylab="£ '000's",nlevels=16,ylim=c(0,20),xlim=c(88,95.5))
lines(seq(88.125,95.5,0.25),average.price,lwd=3)
dev.off()


rm(list = ls(all = TRUE))


Changes<-list()
for (K in c(36,171,353)){  							# Loop over product definitions

	load(file = paste("results",K,".RData",sep=""))
	T <- dim(P)[2]
	J <- dim(Z)[1]
	Pi <- matrix(NA,J,T)
	for (t in 1:T){
		Pi[,t] <- apply(Pi.JK[[t]],1,mean,na.rm=TRUE) # Expected value over draws
	}	
	changers <- changers <- which(apply(1*is.na(Pi),1,sum)!=0)
	Changes[[K]] <- 1*(!is.na(Pi[changers,]))
}

changes <- c(88.875, # 1988q4: adjustable seats and adjustable steering introduced 
			 90.125, # 1990q1: electric sun roofs introduced 
			 90.375, # 1990q2: headlamp cleaners withdrawn;  
			 91.125, # 1991q1: trip computers introduced 
			 92.125, # 1992q1: trip computers withdrawn; 
			 92.625, # 1992q3: trip computers re-introduced for one quarter; 
			 92.875, # 1992q4 radio-only and trip computers withdrawn 
			 93.375, # 1993q2: headlamp cleaners re-introduced 
			 93.875, # 1993q4: air conditioning, power seats and driver’s airbags introduced  
			 94.375, # 1994q2 CD players available 
			 94.875) # 1994q4: headlamp cleaners finally withdrawn 

quarter<-seq(88.125,95.5,0.25)


Product.Indices<-list(list())
Hedonic.Indices<-list(list())
Repack.Indices<-list(list())

for (K in c(36,171,353)){  							# Loop over product definitions

	load(file = paste("results",K,".RData",sep=""))
	
	T <- dim(P)[2]
	J <- dim(Z)[1]
	if (K!=dim(Q)[1]){stop(print("There's a problem with the data pick-up"),call.=TRUE)}
	
	Q[which(is.na(Q))] <- 0 
	Z[which(is.na(Z))] <- 0
	
	if (K==36){
		sigma.characteristics <- 1.611709
		var.sigma.characteristics <- 0.04119106
		sigma.products 		  <- 2.035799 
	} 
	if (K==171){
		sigma.characteristics <- 1.290819
		var.sigma.characteristics <- 0.0181612
		sigma.products 		  <- 1.773457 		
	} 
	if (K==353){
		sigma.characteristics <- 1.500377
		var.sigma.characteristics <- 0.002091437
		sigma.products 		  <- 1.793694 	
	} 

	
	Pi <- matrix(NA,J,T)
	for (t in 1:T){
		Pi[,t] <- apply(Pi.JK[[t]],1,mean,na.rm=TRUE) # Expected value over draws
	}
	#Z[which(Z==0)]<-NA
	colnames(Pi)<-colnames(Z)
	rownames(Pi)<-rownames(Z)
	
	# ------------- #
	# PRODUCT-SPACE #
	# ------------- #
	
	Laspeyres			<- matrix(NA,T,T)
	Paasche 			<- matrix(NA,T,T)
	Fisher 				<- matrix(NA,T,T)
	Sato.Vartia   		<- matrix(NA,T,T)
	Sato.Vartia.Feenstra<- matrix(NA,T,T) 
	for (base in (1:T)){
		for (t in 1:T){
			
			matched.models <- which((Q[,base]>0) & (Q[,t]>0))
			P0 <- P[matched.models,base]		
			P1 <- P[matched.models,t]		
			Q0 <- Q[matched.models,base]		
			Q1 <- Q[matched.models,t]	
			W0 <- (P0 * Q0)/rep(t(P0) %*% Q0,length(matched.models))
			W1 <- (P1 * Q1)/rep(t(P1) %*% Q1,length(matched.models))
			if (all(W0==W1)){
				w <- W0	
			} else {
				w  <- (W1-W0)/(log(W1)-log(W0))
			}
			w <- w/sum(w)
			
			h.1 <- sum(P1*Q1)/sum(P[,t]*Q[,t],na.rm=T)
			h.0 <- sum(P0*Q0)/sum(P[,base]*Q[,base],na.rm=T)
			
			if (length(matched.models)>1){				
				Laspeyres[base,t] <- (t(P1)%*%Q0)/(t(P0)%*%Q0)
				Paasche[base,t] <- (t(P1)%*%Q1)/(t(P0)%*%Q1)
				Fisher[base,t] <- sqrt(Laspeyres[base,t]*Paasche[base,t])
				Sato.Vartia[base,t] <-  prod((P1/P0)^w)	
				Sato.Vartia.Feenstra[base,t] <- ((h.1/h.0)^(1/(sigma.products-1)))*Sato.Vartia[base,t]	
			}
			if (length(matched.models)==1){	# If there is only one matched model		
				Laspeyres[base,t] <- P1/P0
				Paasche[base,t] <- P1/P0
				Fisher[base,t] <- P1/P0
				Sato.Vartia[base,t]   <-  P1/P0	
				Sato.Vartia.Feenstra[base,t] <- ((h.1/h.0)^(1/(sigma.products-1)))*Sato.Vartia[base,t]
			}
		}
	}
	
	sel<-seq(from=T+1,to=T*T-1,length.out=T-1)
	
	Chained.Laspeyres 			<- c(1,cumprod(Laspeyres[sel]))
	Chained.Paasche 			<- c(1,cumprod(Paasche[sel]))
	Chained.Fisher 				<- c(1,cumprod(Fisher[sel]))
	Chained.Sato.Vartia 		<- c(1,cumprod(Sato.Vartia[sel]))
	Chained.Sato.Vartia.Feenstra<- c(1,cumprod(Sato.Vartia.Feenstra[sel]))
	
	Indices <- list(Laspeyres=Laspeyres,Paasche=Paasche,Fisher=Fisher,Sato.Vartia=Sato.Vartia,Sato.Vartia.Feenstra=Sato.Vartia.Feenstra,
	Chained.Laspeyres=Chained.Laspeyres,Chained.Paasche=Chained.Paasche,Chained.Fisher=Chained.Fisher,Chained.Sato.Vartia=Chained.Sato.Vartia,Chained.Sato.Vartia.Feenstra=Chained.Sato.Vartia.Feenstra)
	Product.Indices[[K]]<-Indices
	
	# -------------------- #
	# CHARACTERISTIC-SPACE #
	# -------------------- #
	
	Laspeyres			<- matrix(NA,T,T)
	Paasche 			<- matrix(NA,T,T)
	Fisher 				<- matrix(NA,T,T)
	Sato.Vartia   		<- matrix(NA,T,T)
	Sato.Vartia.Feenstra<- matrix(NA,T,T) 
	se.Laspeyres			<- matrix(NA,T,T)
	se.Paasche 				<- matrix(NA,T,T)
	se.Fisher 				<- matrix(NA,T,T)
	se.Sato.Vartia   		<- matrix(NA,T,T)
	se.Sato.Vartia.Feenstra <- matrix(NA,T,T) 
	
	for (base in (1:T)){
		for (t in 1:T){
			
			matched.characteristics <- which((Z[,base]>0) & (Z[,t]>0))
			Pi0 <- Pi[matched.characteristics,base]		
			Pi1 <- Pi[matched.characteristics,t]		
			Z0 <- Z[matched.characteristics,base]		
			Z1 <- Z[matched.characteristics,t]	
			W0 <- (Pi0 * Z0)/rep(t(Pi0) %*% Z0,length(matched.characteristics))
			W1 <- (Pi1 * Z1)/rep(t(Pi1) %*% Z1,length(matched.characteristics))
			
			if (all(W0==W1)){
				w <- W0	
			} else {
				w  <- (W1-W0)/(log(W1)-log(W0))
				if (TRUE %in% is.na(w)){stop(print("One of the shadow budget shares is fixed",call.=TRUE))}
			}			
			w <- w/sum(w)
			
			eta.1 <- sum(Pi1*Z1)/sum(Pi[,t]*Z[,t],na.rm=T)
			eta.0 <- sum(Pi0*Z0)/sum(Pi[,base]*Z[,base],na.rm=T)
			
			if (length(matched.characteristics)>1){				
				Laspeyres[base,t] <- (t(Pi1)%*%Z0)/(t(Pi0)%*%Z0)
				Paasche[base,t] <- (t(Pi1)%*%Z1)/(t(Pi0)%*%Z1)
				Fisher[base,t] <- sqrt(Laspeyres[base,t]*Paasche[base,t])
				Sato.Vartia[base,t] <-  prod((Pi1/Pi0)^w)	
				Sato.Vartia.Feenstra[base,t] <- ((eta.1/eta.0)^(1/(sigma.characteristics-1)))*Sato.Vartia[base,t]					
				
				V <- as.matrix(bdiag(Var.Pi[[t]][matched.characteristics,matched.characteristics],Var.Pi[[base]][matched.characteristics,matched.characteristics]))
				
				delta <- c(Z0/sum(Pi0*Z0),-sum(Pi0*Z0)*(sum(Pi0*Z0)^-2)*Z0)
				se.Laspeyres[base,t] <- sqrt(delta%*%V%*%delta)
				
				delta <- c(Z1/sum(Pi0*Z1),-sum(Pi0*Z1)*(sum(Pi0*Z1)^-2)*Z1)
				se.Paasche[base,t] <- sqrt(delta%*%V%*%delta)
				
				delta <- c( (0.5*((sum(Pi1*Z0)/sum(Pi0*Z0))*(sum(Pi1*Z1)/sum(Pi0*Z1)))^(-0.5))*(Z0/sum(Pi0*Z1))*(Z1/sum(Pi0*Z1)) , 
						    (0.5*((sum(Pi1*Z0)/sum(Pi0*Z0))*(sum(Pi1*Z1)/sum(Pi0*Z1)))^(-0.5))*(-sum(Pi0*Z0)*(sum(Pi0*Z0)^-2)*Z0)*(-sum(Pi0*Z0)*(sum(Pi0*Z0)^-2)*Z0) )					    
				se.Fisher[base,t] <- sqrt(delta%*%V%*%delta)
				
				delta. <- c(Sato.Vartia[base,t]*(-w*Pi0),Sato.Vartia[base,t]*(w*Pi1))
				se.Sato.Vartia[base,t] <- sqrt(delta%*%V%*%delta)

				delta <- c(-(1/((sigma.characteristics-1)^2))*log(eta.1/eta.0), (eta.1/eta.0)^(1/(sigma.characteristics-1)) )
				V <- diag(c(var.sigma.characteristics,se.Sato.Vartia[base,t]^2))
				se.Sato.Vartia.Feenstra[base,t] <- sqrt(delta%*%V%*%delta)
			
			}
			if (length(matched.models)==1){			
				Laspeyres[base,t] <- Pi1/Pi0
				Paasche[base,t] <- Pi1/Pi0
				Fisher[base,t] <- Pi1/Pi0
				Sato.Vartia[base,t]   <-  Pi1/Pi0	
				Sato.Vartia.Feenstra[base,t] <- ((eta.1/eta.0)^(1/(sigma.characteristics-1)))*Sato.Vartia[base,t]
			}
		}
	}
	diag(se.Laspeyres)<-0
	diag(se.Paasche)<-0
	diag(se.Fisher)<-0
	diag(se.Sato.Vartia)<-0
	diag(se.Sato.Vartia.Feenstra)<-0

	sel<-seq(from=T+1,to=T*T-1,length.out=T-1)
	Chained.Laspeyres 	     	<- c(1,cumprod(Laspeyres[sel]))
	Chained.Paasche 			<- c(1,cumprod(Paasche[sel]))
	Chained.Fisher 				<- c(1,cumprod(Fisher[sel]))
	Chained.Sato.Vartia 		<- c(1,cumprod(Sato.Vartia[sel]))
	Chained.Sato.Vartia.Feenstra<- c(1,cumprod(Sato.Vartia.Feenstra[sel]))
	
	V <- diag(c(0,se.Laspeyres[sel]))
	delta <- Chained.Laspeyres/c(1,Laspeyres[sel])
	se.Chained.Laspeyres <- diag(sqrt(delta*V*delta))
	se.Chained.Laspeyres <- c(0,se.Laspeyres[sel])

	V <- diag(c(0,se.Paasche[sel]))
	delta <- Chained.Paasche/c(1,Paasche[sel])
	se.Chained.Paasche <- diag(sqrt(delta*V*delta))
	se.Chained.Paasche <- c(0,se.Paasche[sel])

	V <- diag(c(0,se.Fisher[sel]))
	delta <- Chained.Fisher/c(1,Fisher[sel])
	se.Chained.Fisher <- diag(sqrt(delta*V*delta))
	se.Chained.Fisher <- c(0,se.Fisher[sel])
	
	V <- diag(c(0,se.Sato.Vartia[sel]))
	delta <- Chained.Sato.Vartia/c(1,Sato.Vartia[sel])
	se.Chained.Sato.Vartia <- diag(sqrt(delta*V*delta))
	se.Chained.Sato.Vartia <- c(0,se.Sato.Vartia[sel])

	
	V <- diag(c(0,se.Sato.Vartia.Feenstra[sel]))
	delta <- Chained.Sato.Vartia.Feenstra/c(1,Sato.Vartia.Feenstra[sel])
	se.Chained.Sato.Vartia.Feenstra <- c(0,diag(sqrt(delta*V*delta)))
	se.Chained.Sato.Vartia.Feenstra <- c(0,se.Sato.Vartia.Feenstra[sel])	
		
	Indices <- list(Laspeyres=Laspeyres,Paasche=Paasche,Fisher=Fisher,Sato.Vartia=Sato.Vartia,Sato.Vartia.Feenstra=Sato.Vartia.Feenstra,
	Chained.Laspeyres=Chained.Laspeyres,Chained.Paasche=Chained.Paasche,Chained.Fisher=Chained.Fisher,Chained.Sato.Vartia=Chained.Sato.Vartia,Chained.Sato.Vartia.Feenstra=Chained.Sato.Vartia.Feenstra,
	se.Laspeyres=se.Laspeyres,se.Paasche=se.Paasche,se.Fisher=se.Fisher,se.Sato.Vartia=se.Sato.Vartia,se.Sato.Vartia.Feenstra=se.Sato.Vartia.Feenstra,
	se.Chained.Laspeyres=se.Chained.Laspeyres,se.Chained.Paasche=se.Chained.Paasche,se.Chained.Fisher=se.Chained.Fisher,
	se.Chained.Sato.Vartia=se.Chained.Sato.Vartia,se.Chained.Sato.Vartia.Feenstra=se.Chained.Sato.Vartia.Feenstra)	
	Hedonic.Indices[[K]] <- Indices	
	
	# Repackaging Model
	
	Repack.Indices[[K]] <- exp(Repack.Index-Repack.Index[,1])

}

# --------------------- #
# Chained Index Numbers #
# --------------------- #
cols <- c("black","purple","blue","red","darkgreen")		
cols <- c("grey50", "grey50", "black", "black",  "grey80")
sym  <- c(3,5,16,8,4)
quartz(width=9,height=7)
pdf(file = paste("..//Manuscript//Figures//Figure_6.pdf",sep=""), onefile = TRUE,width=9,height=7)												# FIGURE 6 #
	plot(NA,NA,typ="l",ylim=c(0.5,1.6),xlim = c(88,95.5), xlab="",ylab="Index")
	for (t in 1:length(changes)){abline(v=changes[t],col="grey71",lty=3)}
	polygon.Product.Sato.Vartia <- NULL
	polygon.Product.Sato.Vartia.Feenstra <- NULL	
	polygon.Hedonic.Sato.Vartia <- NULL
	polygon.Hedonic.Sato.Vartia.Feenstra <- NULL
	polygon.Repackaging <- NULL
	for (K in c(36,171,353)){
		if (K==36){this.line.type = 1}
		if (K==171){this.line.type = 3}
		if (K==353){this.line.type = 2}

		lines(quarter,Product.Indices[[K]]$Chained.Sato.Vartia.Feenstra,col=cols[1],lty=this.line.type)
	   points(quarter,Product.Indices[[K]]$Chained.Sato.Vartia.Feenstra,col=cols[1],pch=sym[1],cex=0.5)
			
		lines(quarter,Product.Indices[[K]]$Chained.Sato.Vartia,col=cols[2],lty=this.line.type)
	   points(quarter,Product.Indices[[K]]$Chained.Sato.Vartia,col=cols[2],pch=sym[2],cex=0.5)
		
		lines(quarter,Hedonic.Indices[[K]]$Chained.Sato.Vartia.Feenstra,col=cols[3],lty=this.line.type)
	   points(quarter,Hedonic.Indices[[K]]$Chained.Sato.Vartia.Feenstra,col=cols[3],pch=sym[3],cex=0.5)	
			
		lines(quarter,Hedonic.Indices[[K]]$Chained.Sato.Vartia,col=cols[4],lty=this.line.type)
	   points(quarter,Hedonic.Indices[[K]]$Chained.Sato.Vartia,col=cols[4],pch=sym[4],cex=0.5)	

		lines(quarter,Repack.Indices[[K]][30,],lty=this.line.type,col=cols[5])
       points(quarter,Repack.Indices[[K]][30,],col=cols[5],pch=sym[5],cex=0.5)	
		
		
		polygon.Product.Sato.Vartia <- rbind(polygon.Product.Sato.Vartia,Product.Indices[[K]]$Chained.Sato.Vartia)

		x.out = sort(c(quarter,seq(min(quarter+0.001),max(quarter)-0.001,length.out=100)))
		
		polygon.Product.Sato.Vartia.Feenstra <- rbind(polygon.Product.Sato.Vartia.Feenstra,
		approx(x=quarter,
		       y=Product.Indices[[K]]$Chained.Sato.Vartia.Feenstra,x.out,
		       method="linear")$y)
		
		polygon.Hedonic.Sato.Vartia <- rbind(polygon.Hedonic.Sato.Vartia,Hedonic.Indices[[K]]$Chained.Sato.Vartia)
		polygon.Hedonic.Sato.Vartia.Feenstra <- rbind(polygon.Hedonic.Sato.Vartia.Feenstra,Hedonic.Indices[[K]]$Chained.Sato.Vartia.Feenstra)
		polygon.Repackaging <- rbind(polygon.Repackaging,Repack.Indices[[K]][30,])	
	}
	
	legend("topleft",c(   "Product-based, mix-adjusted",
						  "Product-based, matched-model",
						  "Hedonic, mix-adjusted",
						  "Hedonic, matched-characteristic",
						  "Time dummy",
						  rep("",17),
						  "Broad product definition (K=36)",
						  "Intermediate product definition (K=171)",
						  "Narrow product definition (K=353)"),
						  #text.col = c("black","purple","blue","red","darkgreen",rep(NA,17),rep("black",3)),
						  text.col = c(cols,rep(NA,17),rep("black",3)),
						  col = c(cols,rep(NA,17),rep("black",3)),
						  pch = c(sym,rep(NA,20)),
						  pt.cex=0.5,
						  lty=c(rep(NA,5),rep(NA,17),1,3,4),
						  bty="n",ncol=1)
	


	poly.y <- c(apply(polygon.Product.Sato.Vartia.Feenstra,2,min,na.rm=TRUE),rev(apply(polygon.Product.Sato.Vartia.Feenstra,2,max,na.rm=TRUE)))
	poly.x <- c(sort(c(quarter,seq(min(quarter+0.001),max(quarter)-0.001,length.out=100))),rev(sort(c(quarter,seq(min(quarter+0.001),max(quarter)-0.001,length.out=100)))))
	polygon(poly.x,poly.y,
	col=rgb(col2rgb(cols[1])[1],col2rgb(cols[1])[2],col2rgb(cols[1])[3],max=255,alpha=(100-90)*255/100),border=NA)

	poly.y <- c(apply(polygon.Product.Sato.Vartia,2,min,na.rm=TRUE),rev(apply(polygon.Product.Sato.Vartia,2,max,na.rm=TRUE)))
	poly.x <- c(quarter,rev(quarter))
	polygon(poly.x,poly.y,
	col=rgb(col2rgb(cols[2])[1],col2rgb(cols[2])[2],col2rgb(cols[2])[3],max=255,alpha=(100-90)*255/100),border=NA)

	poly.y <- c(apply(polygon.Hedonic.Sato.Vartia.Feenstra,2,min,na.rm=TRUE),rev(apply(polygon.Hedonic.Sato.Vartia.Feenstra,2,max,na.rm=TRUE)))
	poly.x <- c(quarter,rev(quarter))
	polygon(poly.x,poly.y,
	col=rgb(col2rgb(cols[3])[1],col2rgb(cols[3])[2],col2rgb(cols[3])[3],max=255,alpha=(100-90)*255/100),border=NA)

	poly.y <- c(apply(polygon.Hedonic.Sato.Vartia,2,min,na.rm=TRUE),rev(apply(polygon.Hedonic.Sato.Vartia,2,max,na.rm=TRUE)))
	poly.x <- c(quarter,rev(quarter))
	polygon(poly.x,poly.y,
	col=rgb(col2rgb(cols[4])[1],col2rgb(cols[4])[2],col2rgb(cols[4])[3],max=255,alpha=(100-90)*255/100),border=NA)
		
	poly.y <- c(apply(polygon.Repackaging,2,min,na.rm=TRUE),rev(apply(polygon.Repackaging,2,max,na.rm=TRUE)))
	poly.x <- c(quarter,rev(quarter))
	polygon(poly.x,poly.y,
	col=rgb(col2rgb(cols[5])[1],col2rgb(cols[5])[2],col2rgb(cols[5])[3],max=255,alpha=(100-90)*255/100),border=NA)	
	
dev.off()
graphics.off()

# ------------------------ #
# Fixed base Index Numbers #
# ------------------------ #
base <- 1
quarter<-seq(88.125,95.5,0.25)
quartz(width=9,height=7)
pdf(file = paste("..//Manuscript//Figures//Figure_3.pdf",sep=""), onefile = TRUE,width=9,height=7)												# FIGURE 3 #
	plot(NA,NA,typ="l",ylim=c(0.5,1.6),xlim = c(88,95.5),xlab="", ylab="Index")
	for (t in 1:length(changes)){abline(v=changes[t],col="grey71",lty=3)}
	polygon.Product.Sato.Vartia <- NULL
	polygon.Product.Sato.Vartia.Feenstra <- NULL	
	polygon.Hedonic.Sato.Vartia <- NULL
	polygon.Hedonic.Sato.Vartia.Feenstra <- NULL
	polygon.Repackaging <- NULL
	for (K in c(36,171,353)){
		if (K==36){this.line.type = 1}
		if (K==171){this.line.type = 3}
		if (K==353){this.line.type = 4}

	lines(quarter,Product.Indices[[K]]$Sato.Vartia.Feenstra[base,],col=cols[1],lty=this.line.type)
	   points(quarter,Product.Indices[[K]]$Sato.Vartia.Feenstra[base,],col=cols[1],pch=sym[1],cex=0.5)
			
		lines(quarter,Product.Indices[[K]]$Sato.Vartia[base,],col=cols[2],lty=this.line.type)
	   points(quarter,Product.Indices[[K]]$Sato.Vartia[base,],col=cols[2],,pch=sym[2],cex=0.5)
		
		lines(quarter,Hedonic.Indices[[K]]$Sato.Vartia.Feenstra[base,],col=cols[3],lty=this.line.type)
	   points(quarter,Hedonic.Indices[[K]]$Sato.Vartia.Feenstra[base,],col=cols[3],pch=sym[3],cex=0.5)	
			
		lines(quarter,Hedonic.Indices[[K]]$Sato.Vartia[base,],col=cols[4],lty=this.line.type)
	   points(quarter,Hedonic.Indices[[K]]$Sato.Vartia[base,],col=cols[4],pch=sym[4],cex=0.5)	

		lines(quarter,Repack.Indices[[K]][30,],lty=this.line.type,col=cols[5])
       points(quarter,Repack.Indices[[K]][30,],col=cols[5],pch=sym[5],cex=0.5)	
	

		
		polygon.Product.Sato.Vartia <- rbind(polygon.Product.Sato.Vartia,Product.Indices[[K]]$Sato.Vartia[base,])

		x.out = sort(c(quarter,seq(min(quarter+0.001),max(quarter)-0.001,length.out=100)))
		
		polygon.Product.Sato.Vartia.Feenstra <- rbind(polygon.Product.Sato.Vartia.Feenstra,
		approx(x=quarter,
		       y=Product.Indices[[K]]$Sato.Vartia.Feenstra[base,],x.out,
		       method="linear")$y)
		
		
		polygon.Hedonic.Sato.Vartia <- rbind(polygon.Hedonic.Sato.Vartia,Hedonic.Indices[[K]]$Sato.Vartia[base,])
		polygon.Hedonic.Sato.Vartia.Feenstra <- rbind(polygon.Hedonic.Sato.Vartia.Feenstra,Hedonic.Indices[[K]]$Sato.Vartia.Feenstra[base,])
		polygon.Repackaging <- rbind(polygon.Repackaging,Repack.Indices[[K]][30,])		
		
	}
						  
		legend("topleft",c(   "Product-based, mix-adjusted",
						  "Product-based, matched-model",
						  "Hedonic, mix-adjusted",
						  "Hedonic, matched-characteristic",
						  "Time dummy",
						  rep("",17),
						  "Broad product definition (K=36)",
						  "Intermediate product definition (K=171)",
						  "Narrow product definition (K=353)"),
						  text.col = c(cols,rep(NA,17),rep("black",3)),
						  col = c(cols,rep(NA,17),rep("black",3)),
						  pch = c(sym,rep(NA,20)),
						  pt.cex=0.5,
						  lty=c(rep(NA,5),rep(NA,17),1,3,4),
						  bty="n",ncol=1)
					  
	poly.y <- c(apply(polygon.Product.Sato.Vartia.Feenstra,2,min,na.rm=TRUE),rev(apply(polygon.Product.Sato.Vartia.Feenstra,2,max,na.rm=TRUE)))
	poly.x <- c(sort(c(quarter,seq(min(quarter+0.001),max(quarter)-0.001,length.out=100))),rev(sort(c(quarter,seq(min(quarter+0.001),max(quarter)-0.001,length.out=100)))))

	polygon(poly.x,poly.y,
	col=rgb(col2rgb(cols[1])[1],col2rgb(cols[1])[2],col2rgb(cols[1])[3],max=255,alpha=(100-90)*255/100),border=NA)

	poly.y <- c(apply(polygon.Product.Sato.Vartia,2,min,na.rm=TRUE),rev(apply(polygon.Product.Sato.Vartia,2,max,na.rm=TRUE)))
	poly.x <- c(quarter,rev(quarter))
	polygon(poly.x,poly.y,
	col=rgb(col2rgb(cols[2])[1],col2rgb(cols[2])[2],col2rgb(cols[2])[3],max=255,alpha=(100-90)*255/100),border=NA)

	poly.y <- c(apply(polygon.Hedonic.Sato.Vartia.Feenstra,2,min,na.rm=TRUE),rev(apply(polygon.Hedonic.Sato.Vartia.Feenstra,2,max,na.rm=TRUE)))
	poly.x <- c(quarter,rev(quarter))
	polygon(poly.x,poly.y,
	col=rgb(col2rgb(cols[3])[1],col2rgb(cols[3])[2],col2rgb(cols[3])[3],max=255,alpha=(100-90)*255/100),border=NA)
	
	poly.y <- c(apply(polygon.Hedonic.Sato.Vartia,2,min,na.rm=TRUE),rev(apply(polygon.Hedonic.Sato.Vartia,2,max,na.rm=TRUE)))
	poly.x <- c(quarter,rev(quarter))
	polygon(poly.x,poly.y,
	col=rgb(col2rgb(cols[4])[1],col2rgb(cols[4])[2],col2rgb(cols[4])[3],max=255,alpha=(100-90)*255/100),border=NA)
	
	poly.y <- c(apply(polygon.Repackaging,2,min,na.rm=TRUE),rev(apply(polygon.Repackaging,2,max,na.rm=TRUE)))
	poly.x <- c(quarter,rev(quarter))
	polygon(poly.x,poly.y,
	col=rgb(col2rgb(cols[5])[1],col2rgb(cols[5])[2],col2rgb(cols[5])[3],max=255,alpha=(100-90)*255/100),border=NA)	

dev.off()
graphics.off()


#############
# Inflation #
#############


quartz(width=9,height=7)
pdf(file = paste("..//Manuscript//Figures//Figure_7.pdf",sep=""), onefile = TRUE,width=9,height=7)												# FIGURE 7 #
	plot(NA,NA,typ="l",xlab="",ylim=c(-30,20),xlim = c(88.75,95.5), ylab="Annual Inflation (pp)")
	for (t in 1:length(changes)){abline(v=changes[t],col="grey71",lty=3)} 
	polygon.Product.Sato.Vartia <- NULL
	polygon.Hedonic.Sato.Vartia <- NULL
	polygon.Hedonic.Sato.Vartia.Feenstra <- NULL
	polygon.Repackaging <- NULL
	for (K in c(36,171,353)){
		if (K==36){this.line.type = 1}
		if (K==171){this.line.type = 3}
		if (K==353){this.line.type = 4}
		
		#if (K==36){lines(quarter[5:T],(Product.Indices[[K]]$Chained.Sato.Vartia.Feenstra[5:T]/Product.Indices[[K]]$Chained.Sato.Vartia.Feenstra[1:(T-4)] - 1),col="blue",lty=this.line.type)}
		lines(quarter[5:T],(Product.Indices[[K]]$Chained.Sato.Vartia[5:T]/Product.Indices[[K]]$Chained.Sato.Vartia[1:(T-4)] - 1)*100,col=cols[2],lty=this.line.type,pch=sym[2])    		  	
	   points(quarter[5:T],(Product.Indices[[K]]$Chained.Sato.Vartia[5:T]/Product.Indices[[K]]$Chained.Sato.Vartia[1:(T-4)] - 1)*100,col=cols[2],pch=sym[2],cex=0.5)


		lines(quarter[5:T],(Hedonic.Indices[[K]]$Chained.Sato.Vartia.Feenstra[5:T]/Hedonic.Indices[[K]]$Chained.Sato.Vartia.Feenstra[1:(T-4)] - 1)*100,col=cols[3],lty=this.line.type,pch=sym[3])
		points(quarter[5:T],(Hedonic.Indices[[K]]$Chained.Sato.Vartia.Feenstra[5:T]/Hedonic.Indices[[K]]$Chained.Sato.Vartia.Feenstra[1:(T-4)] - 1)*100,col=cols[3],pch=sym[3],cex=0.5)

		lines(quarter[5:T],(Hedonic.Indices[[K]]$Chained.Sato.Vartia[5:T]/Hedonic.Indices[[K]]$Chained.Sato.Vartia[1:(T-4)] - 1)*100,col=cols[4],lty=this.line.type,pch=sym[4])
		points(quarter[5:T],(Hedonic.Indices[[K]]$Chained.Sato.Vartia[5:T]/Hedonic.Indices[[K]]$Chained.Sato.Vartia[1:(T-4)] - 1)*100,col=cols[4],pch=sym[4],cex=0.5)

		lines(quarter[5:T],(Repack.Indices[[K]][30,][5:T]/Repack.Indices[[K]][30,][1:(T-4)] - 1)*100,col=cols[5],lty=this.line.type,pch=sym[4])
		points(quarter[5:T],(Repack.Indices[[K]][30,][5:T]/Repack.Indices[[K]][30,][1:(T-4)] - 1)*100,col=cols[5],pch=sym[4],cex=0.5)
		
		polygon.Product.Sato.Vartia <- rbind(polygon.Product.Sato.Vartia,(Product.Indices[[K]]$Chained.Sato.Vartia[5:T]/Product.Indices[[K]]$Chained.Sato.Vartia[1:(T-4)] - 1)*100)
		polygon.Hedonic.Sato.Vartia.Feenstra <- rbind(polygon.Hedonic.Sato.Vartia.Feenstra,(Hedonic.Indices[[K]]$Chained.Sato.Vartia.Feenstra[5:T]/Hedonic.Indices[[K]]$Chained.Sato.Vartia.Feenstra[1:(T-4)] - 1)*100)
		polygon.Hedonic.Sato.Vartia <- rbind(polygon.Hedonic.Sato.Vartia,(Hedonic.Indices[[K]]$Chained.Sato.Vartia[5:T]/Hedonic.Indices[[K]]$Chained.Sato.Vartia[1:(T-4)] - 1)*100)
		polygon.Repackaging <- rbind(polygon.Repackaging,(Repack.Indices[[K]][30,][5:T]/Repack.Indices[[K]][30,][1:(T-4)] - 1)*100)	

	}
	
	poly.x <- c(quarter[5:T],rev(quarter[5:T]))
	
	poly.y <- c(apply(polygon.Product.Sato.Vartia,2,min,na.rm=TRUE),rev(apply(polygon.Product.Sato.Vartia,2,max,na.rm=TRUE)))
	polygon(poly.x,poly.y,
	col=rgb(col2rgb(cols[2])[1],col2rgb(cols[2])[2],col2rgb(cols[2])[3],max=255,alpha=(100-90)*255/100),border=NA)

	poly.y <- c(apply(polygon.Hedonic.Sato.Vartia.Feenstra,2,min,na.rm=TRUE),rev(apply(polygon.Hedonic.Sato.Vartia.Feenstra,2,max,na.rm=TRUE)))
	polygon(poly.x,poly.y,
	col=rgb(col2rgb(cols[3])[1],col2rgb(cols[3])[2],col2rgb(cols[3])[3],max=255,alpha=(100-90)*255/100),border=NA)	

	poly.y <- c(apply(polygon.Hedonic.Sato.Vartia,2,min,na.rm=TRUE),rev(apply(polygon.Hedonic.Sato.Vartia,2,max,na.rm=TRUE)))
	polygon(poly.x,poly.y,
	col=rgb(col2rgb(cols[4])[1],col2rgb(cols[4])[2],col2rgb(cols[4])[3],max=255,alpha=(100-90)*255/100),border=NA)
	
	poly.y <- c(apply(polygon.Repackaging,2,min,na.rm=TRUE),rev(apply(polygon.Repackaging,2,max,na.rm=TRUE)))
	polygon(poly.x,poly.y,
	col=rgb(col2rgb(cols[5])[1],col2rgb(cols[5])[2],col2rgb(cols[5])[3],max=255,alpha=(100-90)*255/100),border=NA)	
	
	legend("topleft",c("Product-based, matched-model",
					  #"Product-based, mix-adjusted",
					  "Hedonic, mix-adjusted",
					  "Hedonic, matched-characteristic",
					  "Time dummy",
					  rep("",17),
					  "Broad product definition (K=36)",
					  "Intermediate product definition (K=171)",
					  "Narrow product definition (K=353)"),
					  text.col = c(cols[2:5],rep(NA,17),rep("black",3)),
					  col = c(cols[2:5],rep(NA,17),rep("black",3)),
					  pch = c(sym[2:5],rep(NA,20)),
  					  pt.cex=0.5,
 					  lty=c(rep(NA,4),rep(NA,17),1,3,4),
					  bty="n",ncol=1)					
	
dev.off()
graphics.off()


quartz(width=9,height=7)
pdf(file = paste("..//Manuscript//Figures//Figure_4.pdf",sep=""), onefile = TRUE,width=9,height=7)   											# FIGURE 4 #
	plot(NA,NA,typ="l",xlab="",ylim=c(-30,20),xlim = c(88.75,95.5), ylab="Annual Inflation (pp)")
	for (t in 1:length(changes)){abline(v=changes[t],col="grey71",lty=3)}
	polygon.Product.Sato.Vartia <- NULL
	polygon.Hedonic.Sato.Vartia <- NULL
	polygon.Hedonic.Sato.Vartia.Feenstra <- NULL
	polygon.Repackaging <- NULL
	
	for (K in c(36,171,353)){
		if (K==36){this.line.type = 1}
		if (K==171){this.line.type = 3}
		if (K==353){this.line.type = 2}
		lines(quarter[5:T],(Product.Indices[[K]]$Sato.Vartia[base,5:T]/Product.Indices[[K]]$Sato.Vartia[base,1:(T-4)] - 1)*100,col=cols[2],lty=this.line.type)
		points(quarter[5:T],(Product.Indices[[K]]$Sato.Vartia[base,5:T]/Product.Indices[[K]]$Sato.Vartia[base,1:(T-4)] - 1)*100,col=cols[2],pch=sym[2],cex=0.5)
		
		lines(quarter[5:T],(Hedonic.Indices[[K]]$Sato.Vartia.Feenstra[base,5:T]/Hedonic.Indices[[K]]$Sato.Vartia.Feenstra[base,1:(T-4)] - 1)*100,col=cols[3],lty=this.line.type)
		points(quarter[5:T],(Hedonic.Indices[[K]]$Sato.Vartia.Feenstra[base,5:T]/Hedonic.Indices[[K]]$Sato.Vartia.Feenstra[base,1:(T-4)] - 1)*100,col=cols[3],pch=sym[3],cex=0.5)

		lines(quarter[5:T],(Hedonic.Indices[[K]]$Sato.Vartia[base,5:T]/Hedonic.Indices[[K]]$Sato.Vartia[base,1:(T-4)] - 1)*100,col=cols[4],lty=this.line.type)
		points(quarter[5:T],(Hedonic.Indices[[K]]$Sato.Vartia[base,5:T]/Hedonic.Indices[[K]]$Sato.Vartia[base,1:(T-4)] - 1)*100,col=cols[4],pch=sym[4],cex=0.5)
		
		lines(quarter[5:T],(Repack.Indices[[K]][30,][5:T]/Repack.Indices[[K]][30,][1:(T-4)] - 1)*100,col=cols[5],lty=this.line.type)
		points(quarter[5:T],(Repack.Indices[[K]][30,][5:T]/Repack.Indices[[K]][30,][1:(T-4)] - 1)*100,col=cols[5],pch=sym[5],cex=0.5)

		
		polygon.Product.Sato.Vartia <- rbind(polygon.Product.Sato.Vartia,(Product.Indices[[K]]$Sato.Vartia[base,5:T]/Product.Indices[[K]]$Sato.Vartia[base,1:(T-4)] - 1)*100)
		polygon.Hedonic.Sato.Vartia <- rbind(polygon.Hedonic.Sato.Vartia,(Hedonic.Indices[[K]]$Sato.Vartia[base,5:T]/Hedonic.Indices[[K]]$Sato.Vartia[base,1:(T-4)] - 1)*100)
		polygon.Hedonic.Sato.Vartia.Feenstra <- rbind(polygon.Hedonic.Sato.Vartia.Feenstra,(Hedonic.Indices[[K]]$Sato.Vartia.Feenstra[base,5:T]/Hedonic.Indices[[K]]$Sato.Vartia.Feenstra[base,1:(T-4)] - 1)*100)
		polygon.Repackaging <- rbind(polygon.Repackaging,(Repack.Indices[[K]][30,][5:T]/Repack.Indices[[K]][30,][1:(T-4)] - 1)*100)	
	}

	poly.x <- c(quarter[5:T],rev(quarter[5:T]))
	
	poly.y <- c(apply(polygon.Product.Sato.Vartia,2,min,na.rm=TRUE),rev(apply(polygon.Product.Sato.Vartia,2,max,na.rm=TRUE)))
	polygon(poly.x,poly.y,
	col=rgb(col2rgb(cols[2])[1],col2rgb(cols[2])[2],col2rgb(cols[2])[3],max=255,alpha=(100-90)*255/100),border=NA)
	
	poly.y <- c(apply(polygon.Hedonic.Sato.Vartia.Feenstra,2,min,na.rm=TRUE),rev(apply(polygon.Hedonic.Sato.Vartia.Feenstra,2,max,na.rm=TRUE)))
	polygon(poly.x,poly.y,
	col=rgb(col2rgb(cols[3])[1],col2rgb(cols[3])[2],col2rgb(cols[3])[3],max=255,alpha=(100-90)*255/100),border=NA)

	poly.y <- c(apply(polygon.Hedonic.Sato.Vartia,2,min,na.rm=TRUE),rev(apply(polygon.Hedonic.Sato.Vartia,2,max,na.rm=TRUE)))
	polygon(poly.x,poly.y,
	col=rgb(col2rgb(cols[4])[1],col2rgb(cols[4])[2],col2rgb(cols[4])[3],max=255,alpha=(100-90)*255/100),border=NA)
	
	poly.y <- c(apply(polygon.Repackaging,2,min,na.rm=TRUE),rev(apply(polygon.Repackaging,2,max,na.rm=TRUE)))
	polygon(poly.x,poly.y,
	col=rgb(col2rgb(cols[5])[1],col2rgb(cols[5])[2],col2rgb(cols[5])[3],max=255,alpha=(100-90)*255/100),border=NA)	
	

	legend("topleft",c("Product-based, matched-model",
					  #"Product-based, mix-adjusted",
					  "Hedonic, mix-adjusted",
					  "Hedonic, matched-characteristic",
					  "Time dummy",
					  rep("",17),
					  "Broad product definition (K=36)",
					  "Intermediate product definition (K=171)",
					  "Narrow product definition (K=353)"),
					  text.col = c(cols[2:5],rep(NA,17),rep("black",3)),
					  col = c(cols[2:5],rep(NA,17),rep("black",3)),
					  pch = c(sym[2:5],rep(NA,20)),
  					  pt.cex=0.5,
 					  lty=c(rep(NA,4),rep(NA,17),1,3,4),
					  bty="n",ncol=1)					

dev.off()	
graphics.off()

# ---------------------------------------------------- #
# Information on the prevalence/value of matched pairs #
# ---------------------------------------------------- #

proportion.matched.models <- list()
value.matched.models <- list()
proportion.matched.characteristics <- list()
value.matched.characteristics <- list()
for (K in c(36,171,353)){
	load(file = paste("results",K,".RData",sep=""))

	T <- dim(P)[2]
	J <- dim(Z)[1]
	
	Q[which(is.na(Q))] <- 0 
	Z[which(is.na(Z))] <- 0	
	
	# Matched models #
	
	prop.matched.models <- matrix(NA,T,T)
	rownames(prop.matched.models)<-paste("Base period = ",1:T)
	colnames(prop.matched.models)<-paste("Period = ",1:T)
	share.matched.models <- prop.matched.models
	for (base in 1:T){
		for (t in 1:T){
			matched.models <- which((Q[,base]>0) & (Q[,t]>0))
			prop.matched.models[base,t] <- length(matched.models)/sum(Q[,base]>0)  # N matched models/N models in the base period.
			share.matched.models[base,t] <- sum(P[matched.models,base]*Q[matched.models,base],na.rm=TRUE)/sum(P[,base]*Q[,base],na.rm=TRUE)
		}
	}
	proportion.matched.models[[K]]<-prop.matched.models
	value.matched.models[[K]]<-share.matched.models

	# Matched characteristics #

	prop.matched.characteristics <- matrix(NA,T,T)	
	rownames(prop.matched.characteristics)<-paste("Base period = ",1:T)
	colnames(prop.matched.characteristics)<-paste("Period = ",1:T)
	share.matched.characteristics <- prop.matched.characteristics
	for (base in 1:T){
		for (t in 1:T){
			matched.characteristics <- which((Z[,base]>0) & (Z[,t]>0))
			prop.matched.characteristics[base,t] <- length(matched.characteristics)/sum(Z[,base]>0) # N matched char's/N char's available in the base period.
			share.matched.characteristics[base,t] <- sum(Pi[matched.characteristics,base]*Z[matched.characteristics,base],na.rm=TRUE)/sum(Pi[,base]*Z[,base],na.rm=TRUE)
		}
	}
	proportion.matched.characteristics[[K]]<-prop.matched.characteristics
	value.matched.characteristics[[K]]<-share.matched.characteristics
	
}


quartz(width=9,height=7)
pdf(file = paste("..//Manuscript//Figures//Figure_2a.pdf",sep=""), onefile = TRUE,width=9,height=7)												# FIGURE 2(a) #
plot(NA,NA,typ="l",xlab="Quarter lag length",ylim=c(0,1),xlim = c(1,T), ylab="Proportion, matched models")
dot.colours <- rep("black",353)
c("grey","pink","cornflowerblue")
line.colours <- c("black","red","blue")
for (K in c(36,171,353)){
	y <- NULL
	x <- NULL
	if (K==36){
		dot.colour = "black"
		line.colour <- "black"
		dot.cex <- 0.5
		dot <- 0
	}
	if (K==171){
		dot.colour = "grey80"
		line.colour <- "grey80"
		dot.cex <- 0.5
		dot<-1
	}
	if (K==353){
		dot.colour = "grey50"
		line.colour <- "grey50"
		dot.cex <- 0.5
		dot<-2
	}	
	
	for (base in 1:T){
			points(abs(1:T-base),proportion.matched.models[[K]][base,],col=dot.colour,pch=dot,cex=dot.cex)
			x <- c(x,abs(1:T-base))
			y <- c(y,proportion.matched.models[[K]][base,])
	}
	
	mean <- NULL
	for (diff in 0:29){mean <- c(mean,mean(y[which(x==diff)]))}
	lines(0:29,mean,col=line.colour)
}

dev.off()

quartz(width=9,height=7)
pdf(file = paste("..//Manuscript//Figures//Figure_2b.pdf",sep=""), onefile = TRUE,width=9,height=7)                  							# FIGURE 2(b) #
plot(NA,NA,typ="l",xlab="Quarter lag length",ylim=c(0,1),xlim = c(1,T), ylab="Proportion, matched characteristics")
dot.colours <- rep("black",353)
c("grey","pink","cornflowerblue")
line.colours <- c("black","red","blue")
for (K in c(36,171,353)){
	y <- NULL
	x <- NULL
	if (K==36){
		dot.colour = "black"
		line.colour <- "black"
		dot.cex <- 0.66
	}
	if (K==171){
		dot.colour = "pink"
		line.colour <- "red"
		dot.cex <- 0.5
	}
	if (K==353){
		dot.colour = "cornflowerblue"
		line.colour <- "blue"
		dot.cex <- 0.33
	}	
	
	for (base in 1:T){
			points(abs(1:T-base),proportion.matched.characteristics[[K]][base,],col="black",pch=16,cex=0.5)
			x <- c(x,abs(1:T-base))
			y <- c(y,proportion.matched.characteristics[[K]][base,])
	}
	
	mean <- NULL
	for (diff in 0:29){mean <- c(mean,mean(y[which(x==diff)]))}
	lines(0:29,mean,col="black")
}
dev.off()

#----------------------------------------------------------#
# Broad definition, all indices, plus confidence intervals #
#----------------------------------------------------------#

K <- 36

#------------#
# Fixed Base #
#------------#

base <- 1
quarter<-seq(88.125,95.5,0.25)
cols <- c(rep("grey70",4),rep("black",4),"grey20")
quartz(width=9,height=7)
pdf(file = paste("..//Manuscript//Figures//Figure_8.pdf",sep=""), onefile = TRUE,width=9,height=7)                  							# FIGURE 8 #
	plot(NA,NA,typ="l",xlab="Quarter",ylim=c(0.95,1.5),xlim = c(88,95.5), ylab="Index")
	for (t in 1:length(changes)){abline(v=changes[t],col="grey71",lty=3)}
	
	this.line.type = 1
	
	lines(quarter,Product.Indices[[K]]$Laspeyres[base,],col=cols[1],lty=4)
	lines(quarter,Product.Indices[[K]]$Paasche[base,],col=cols[2],lty=3)
	lines(quarter,Product.Indices[[K]]$Fisher[base,],col=cols[3],lty=2)
	lines(quarter,Product.Indices[[K]]$Sato.Vartia[base,],col=cols[4],lty=1)
	
	#lines(quarter,Product.Indices[[K]]$Sato.Vartia.Feenstra[base,],col="purple",lty=this.line.type)
	
	lines(quarter,Hedonic.Indices[[K]]$Laspeyres[base,],col=cols[5],lty=4)
	lines(quarter,Hedonic.Indices[[K]]$Paasche[base,],col=cols[6],lty=3)
	lines(quarter,Hedonic.Indices[[K]]$Fisher[base,],col=cols[7],lty=2)
	lines(quarter,Hedonic.Indices[[K]]$Sato.Vartia[base,],col=cols[8],lty=1)	
	
	lines(quarter,Hedonic.Indices[[K]]$Sato.Vartia.Feenstra[base,],col=cols[9],lty=1)
	points(quarter,Hedonic.Indices[[K]]$Sato.Vartia.Feenstra[base,],col=cols[9],pch=16,cex=0.5)
	
	index <- Hedonic.Indices[[K]]$Laspeyres[base,]
	se.index <- Hedonic.Indices[[K]]$se.Laspeyres[base,]
	conf.int <- cbind(rbind(quarter,index+1.96*se.index),rbind(rev(quarter),rev(index-1.96*se.index)))
	polygon(conf.int[1,],conf.int[2,],col=rgb(col2rgb(cols[5])[1],col2rgb(cols[5])[2],col2rgb(cols[5])[3],max=255,alpha=(100-90)*255/100),border=NA)
	
	index <- Hedonic.Indices[[K]]$Paasche[base,]
	se.index <- Hedonic.Indices[[K]]$se.Paasche[base,]
	conf.int <- cbind(rbind(quarter,index+1.96*se.index),rbind(rev(quarter),rev(index-1.96*se.index)))
	polygon(conf.int[1,],conf.int[2,],col=rgb(col2rgb(cols[6])[1],col2rgb(cols[6])[2],col2rgb(cols[6])[3],max=255,alpha=(100-90)*255/100),border=NA)	

	index <- Hedonic.Indices[[K]]$Fisher[base,]
	se.index <- Hedonic.Indices[[K]]$se.Fisher[base,]
	conf.int <- cbind(rbind(quarter,index+1.96*se.index),rbind(rev(quarter),rev(index-1.96*se.index)))
	polygon(conf.int[1,],conf.int[2,],col=rgb(col2rgb(cols[7])[1],col2rgb(cols[7])[2],col2rgb(cols[7])[3],max=255,alpha=(100-90)*255/100),border=NA)	

	index <- Hedonic.Indices[[K]]$Sato.Vartia[base,]
	se.index <- Hedonic.Indices[[K]]$se.Sato.Vartia[base,]
	conf.int <- cbind(rbind(quarter,index+1.96*se.index),rbind(rev(quarter),rev(index-1.96*se.index)))
	polygon(conf.int[1,],conf.int[2,],col=rgb(col2rgb(cols[8])[1],col2rgb(cols[8])[2],col2rgb(cols[8])[3],max=255,alpha=(100-90)*255/100),border=NA)		

	index <- Hedonic.Indices[[K]]$Sato.Vartia.Feenstra[base,]
	se.index <- Hedonic.Indices[[K]]$se.Sato.Vartia.Feenstra[base,]
	conf.int <- cbind(rbind(quarter,index+1.96*se.index),rbind(rev(quarter),rev(index-1.96*se.index)))
	polygon(conf.int[1,],conf.int[2,],col=rgb(col2rgb(cols[9])[1],col2rgb(cols[9])[2],col2rgb(cols[9])[3],max=255,alpha=(100-90)*255/100),border=NA)	
	
	

legend.text <- c("Product-based Laspeyres","Product-based Paasche","Product-based Fisher","Product-based Sato-Vartia",
				"Hedonic Laspeyres","Hedonic Paasche","Hedonic Fisher","Hedonic Sato-Vartia","Hedonic mix-adjusted Sato-Vartia-Feenstra")
legend.ltype <- c(4:1,4:1,1)
legend("topleft",legend.text,col=cols,lty=legend.ltype,pch = c(rep(NA,8),16),pt.cex=0.5,bty="n")	

#dev.off()
graphics.off()



#---------#	
# Chained #
#---------#
quarter<-seq(88.125,95.5,0.25)
quartz(width=9,height=7)
pdf(file = paste("..//Manuscript//Figures//Figure_9.pdf",sep=""), onefile = TRUE,width=9,height=7)                  							# FIGURE 9 #
	plot(NA,NA,typ="l",xlab="Quarter",ylim=c(0.95,1.5),xlim = c(88,95.5), ylab="Index")
	for (t in 1:length(changes)){abline(v=changes[t],col="grey71",lty=3)}
	
	this.line.type = 1
	
	lines(quarter,Product.Indices[[K]]$Chained.Laspeyres,col=cols[1],lty=4)
	lines(quarter,Product.Indices[[K]]$Chained.Paasche,col=cols[2],lty=3)
	lines(quarter,Product.Indices[[K]]$Chained.Fisher,col=cols[3],lty=2)
	lines(quarter,Product.Indices[[K]]$Chained.Sato.Vartia,col=cols[4],lty=1)
	
	#lines(quarter,Product.Indices[[K]]$Chained.Sato.Vartia.Feenstra,col="purple",lty=this.line.type)
	
	lines(quarter,Hedonic.Indices[[K]]$Chained.Laspeyres,col=cols[5],lty=4)
	lines(quarter,Hedonic.Indices[[K]]$Chained.Paasche,col=cols[6],lty=3)
	lines(quarter,Hedonic.Indices[[K]]$Chained.Fisher,col=cols[7],lty=2)
	lines(quarter,Hedonic.Indices[[K]]$Chained.Sato.Vartia,col=cols[8],lty=1)	
	
	lines(quarter,Hedonic.Indices[[K]]$Chained.Sato.Vartia.Feenstra,col=cols[9],lty=1)
	points(quarter,Hedonic.Indices[[K]]$Chained.Sato.Vartia.Feenstra,col=cols[9],pch=16,cex=0.5)
		

	
	index <- Hedonic.Indices[[K]]$Chained.Laspeyres
	se.index <- Hedonic.Indices[[K]]$se.Chained.Laspeyres
	conf.int <- cbind(rbind(quarter,index+1.96*se.index),rbind(rev(quarter),rev(index-1.96*se.index)))
	polygon(conf.int[1,],conf.int[2,],col=rgb(col2rgb(cols[5])[1],col2rgb(cols[5])[2],col2rgb(cols[5])[3],max=255,alpha=(100-90)*255/100),border=NA)
	
	index <- Hedonic.Indices[[K]]$Chained.Paasche
	se.index <- Hedonic.Indices[[K]]$se.Chained.Paasche
	conf.int <- cbind(rbind(quarter,index+1.96*se.index),rbind(rev(quarter),rev(index-1.96*se.index)))
	polygon(conf.int[1,],conf.int[2,],col=rgb(col2rgb(cols[6])[1],col2rgb(cols[6])[2],col2rgb(cols[6])[3],max=255,alpha=(100-90)*255/100),border=NA)	

	index <- Hedonic.Indices[[K]]$Chained.Fisher
	se.index <- Hedonic.Indices[[K]]$se.Chained.Fisher
	conf.int <- cbind(rbind(quarter,index+1.96*se.index),rbind(rev(quarter),rev(index-1.96*se.index)))
	polygon(conf.int[1,],conf.int[2,],col=rgb(col2rgb(cols[7])[1],col2rgb(cols[7])[2],col2rgb(cols[7])[3],max=255,alpha=(100-90)*255/100),border=NA)	

	index <- Hedonic.Indices[[K]]$Chained.Sato.Vartia
	se.index <- Hedonic.Indices[[K]]$se.Chained.Sato.Vartia
	conf.int <- cbind(rbind(quarter,index+1.96*se.index),rbind(rev(quarter),rev(index-1.96*se.index)))
	polygon(conf.int[1,],conf.int[2,],col=rgb(col2rgb(cols[8])[1],col2rgb(cols[8])[2],col2rgb(cols[8])[3],max=255,alpha=(100-90)*255/100),border=NA)		

	index <- Hedonic.Indices[[K]]$Chained.Sato.Vartia.Feenstra
	se.index <- Hedonic.Indices[[K]]$se.Chained.Sato.Vartia.Feenstra
	conf.int <- cbind(rbind(quarter,index+1.96*se.index),rbind(rev(quarter),rev(index-1.96*se.index)))
	polygon(conf.int[1,],conf.int[2,],col=rgb(col2rgb(cols[9])[1],col2rgb(cols[9])[2],col2rgb(cols[9])[3],max=255,alpha=(100-90)*255/100),border=NA)	
	
	

legend.text <- c("Product-based Laspeyres","Product-based Paasche","Product-based Fisher","Product-based Sato-Vartia",
				"Hedonic Laspeyres","Hedonic Paasche","Hedonic Fisher","Hedonic Sato-Vartia","Hedonic mix-adjusted Sato-Vartia-Feenstra")
legend.cols <- cols
legend.ltype <- c(4:1,4:1,1)
legend("topleft",legend.text,col=cols,lty=legend.ltype,pch = c(rep(NA,8),16),pt.cex=0.5,bty="n")	

dev.off()
graphics.off()

#-----------------#
# Repacking Index #
#-----------------#

# Levels #
quartz(width=9,height=7)
pdf(file = paste("..//Manuscript//Figures//Figure_10.pdf",sep=""), onefile = TRUE,width=9,height=7)                  							# FIGURE 10 #

plot(NA,NA,typ="l",xlab="",ylim=c(0.99,1.35),xlim = c(88,95.5), ylab="Index")	
for (t in 1:length(changes)){abline(v=changes[t],col="grey71",lty=3)}

for (K in c(36,171,353)){
		if (K==36){colour<-"black"}
		if (K==171){colour<-"grey40"}
		if (K==353){colour<-"grey80"}
		if (K==36){dot<-3}
		if (K==171){dot<-15}
		if (K==353){dot<-16}

	for (t in T:1){
		color=rgb(col2rgb(colour)[1],col2rgb(colour)[2],col2rgb(colour)[3],max=255,alpha=round(255*(t/T),0))
		lines(quarter,Repack.Indices[[K]][t,],col=color)
		points(quarter,Repack.Indices[[K]][t,],col=color,pch=dot,cex=0.75)
	}
}
legend.text <- c("Broad definition of products (K=36)","Intermediate definition of products (K=171)","Narrow definition of products (K=353)")
legend.cols <- c("black","grey40","grey80")
legend.ltype <- c(1,1,1)
legend("topleft",legend.text,lty=legend.ltype,pch = c(3,15,16),pt.cex=0.5,bty="n",col=legend.cols)	

dev.off()
graphics.off()

# Inflation #

quartz(width=9,height=7)
pdf(file = paste("..//Manuscript//Figures//Figure_11.pdf",sep=""), onefile = TRUE,width=9,height=7)                  							# FIGURE 11 #
plot(NA,NA,typ="l",xlab="Quarter",ylim=c(-5,10),xlim = c(89,95.5), ylab="Annual Inflation (pp)")	
for (t in 1:length(changes)){abline(v=changes[t],col="grey71",lty=3)}

for (K in c(36,171,353)){
		if (K==36){colour<-"black"}
		if (K==171){colour<-"grey40"}
		if (K==353){colour<-"grey80"}
		if (K==36){dot<-3}
		if (K==171){dot<-15}
		if (K==353){dot<-16}

	for (t in T:1){
		color=rgb(col2rgb(colour)[1],col2rgb(colour)[2],col2rgb(colour)[3],max=255,alpha=round(255*(t/T),0))
		lines(quarter[5:T],(Repack.Indices[[K]][t,][5:T]/Repack.Indices[[K]][t,][1:(T-4)] - 1)*100,col=color)
		points(quarter[5:T],(Repack.Indices[[K]][t,][5:T]/Repack.Indices[[K]][t,][1:(T-4)] - 1)*100,col=color,pch=dot,cex=0.75)
	}	
}
legend.text <- c("Broad definition of products (K=36)","Intermediate definition of products (K=171)","Narrow definition of products (K=353)")
legend.cols <- c("black","grey40","grey80")
legend.ltype <- c(1,1,1)
legend("bottomleft",legend.text,lty=legend.ltype,pch = c(3,15,16),pt.cex=0.5,bty="n",col=legend.cols)	
dev.off()
graphics.off()


# ---------------------#
# Sensitivity to sigma # 
# ---------------------#
sigmas <- seq(from=1.05,to=5,by=0.05)
for (K in c(36,171,353)){
	
	load(file = paste("results",K,".RData",sep=""))
	
	T <- dim(P)[2]
	J <- dim(Z)[1]
	if (K!=dim(Q)[1]){stop(print("There's a problem with the data pick-up"),call.=TRUE)}
	
	Q[which(is.na(Q))] <- 0 
	Z[which(is.na(Z))] <- 0
	
	if (K==36){ sigma.characteristics <- sort(c(1.6117092,sigmas))} 
	if (K==171){sigma.characteristics <- sort(c(1.290819,sigmas))}
	if (K==353){sigma.characteristics <- sort(c(1.500377,sigmas))}

	Chained.Sato.Vartia.Feenstra <- matrix(NA,length(sigma.characteristics),T)
	Fixed.Sato.Vartia.Feenstra <- matrix(NA,length(sigma.characteristics),T)	
	
	Pi <- matrix(NA,J,T)
	for (t in 1:T){
		Pi[,t] <- apply(Pi.JK[[t]],1,mean,na.rm=TRUE) # Expected value over draws
	}
#	Z[which(Z==0)]<-NA
	colnames(Pi)<-colnames(Z)
	rownames(Pi)<-rownames(Z)
	
	for (sig in 1:length(sigma.characteristics)){
	
		for (base in (1:T)){
			for (t in 1:T){
				
				matched.characteristics <- which((Z[,base]>0) & (Z[,t]>0))
				Pi0 <- Pi[matched.characteristics,base]		
				Pi1 <- Pi[matched.characteristics,t]		
				Z0 <- Z[matched.characteristics,base]		
				Z1 <- Z[matched.characteristics,t]	
				W0 <- (Pi0 * Z0)/rep(t(Pi0) %*% Z0,length(matched.characteristics))
				W1 <- (Pi1 * Z0)/rep(t(Pi1) %*% Z1,length(matched.characteristics))

				if (all(W0==W1)){
					w <- W0	
				} else {
					w  <- (W1-W0)/(log(W1)-log(W0))
					if (TRUE %in% is.na(w)){stop(print("One of the budget shares is fixed",call.=TRUE))}
				}
				
				eta.1 <- sum(Pi1*Z1)/sum(Pi[,t]*Z[,t],na.rm=T)
				eta.0 <- sum(Pi0*Z0)/sum(Pi[,base]*Z[,base],na.rm=T)
			
				if (length(matched.characteristics)>1){				
					Sato.Vartia[base,t] <-  prod((Pi1/Pi0)^w)	
					Sato.Vartia.Feenstra[base,t] <- ((eta.1/eta.0)^(1/(sigma.characteristics[sig]-1)))*Sato.Vartia[base,t]	
				}
			}
		}
		
		sel<-seq(from=T+1,to=T*T-1,length.out=T-1)

		Chained.Sato.Vartia 		       <- c(1,cumprod(Sato.Vartia[sel]))
		Chained.Sato.Vartia.Feenstra[sig,] <- c(1,cumprod(Sato.Vartia.Feenstra[sel]))	
		Fixed.Sato.Vartia.Feenstra[sig,]   <- Sato.Vartia.Feenstra[1,]     # Base period = 1 index #
		Fixed.Sato.Vartia                  <- Sato.Vartia[1,]
		
	}
	
	if (K==36){ sigma.characteristics <- sort(c(1.6117092,sigmas))} 
	if (K==171){sigma.characteristics <- sort(c(1.290819,sigmas))}
	if (K==353){sigma.characteristics <- sort(c(1.500377,sigmas))}	
			
	bias <- Chained.Sato.Vartia.Feenstra[,T]/Chained.Sato.Vartia[T]
	if (K==36){
		quartz(width=8,height=5)
		pdf(file = paste("..//Manuscript//Figures//Figure_5.pdf",sep=""), onefile = TRUE,width=9,height=7)                  							# FIGURE 5 #
		
		plot(sigma.characteristics,bias,typ="l",xlab=expression(sigma),ylab="Bias")
		sel <- which(sigma.characteristics==1.6117092)
		points(sigma.characteristics[sel],bias[sel],pch=16)
		text(sigma.characteristics[sel],bias[sel],labels="K=36",pos=1)
		Bias <- bias
	}
	if (K==171){ 
		lines(sigma.characteristics,bias,lty=3)
		sel <- which(sigma.characteristics==1.290819)
		points(sigma.characteristics[sel],bias[sel],pch=16)
		text(sigma.characteristics[sel],bias[sel],labels="K=171",pos=2)
		Bias <- rbind(Bias,bias)
	}
	if (K==353){ 
		lines(sigma.characteristics,bias,lty=4)
		sel <- which(sigma.characteristics==1.500377)
		points(sigma.characteristics[sel],bias[sel],pch=16)
		text(sigma.characteristics[sel],bias[sel],labels="K=353",pos=3)
		Bias <- rbind(Bias,bias)
		
	
	sigma.K <- c("Broad definition (K=36)","Intermediate definition (K=171)","Narrow definition (K=353)")

	legend("bottomright",sigma.K,col=rep("black",3),lty=c(1,3,4),bty="n")
	}
}

dev.off()
graphics.off()

# Some further figures which never made it into the final paper #


quartz(width=9,height=7)
plot(NA,NA,typ="l",xlab="Quarter lag length",ylim=c(0,1),xlim = c(1,T), ylab="Value share, matched models")
dot.colours <- rep("black",353)
c("grey","pink","cornflowerblue")
line.colours <- c("black","red","blue")
for (K in c(36,171,353)){
	y <- NULL
	x <- NULL
	if (K==36){
		dot.colour = "grey"
		line.colour <- "black"
		dot.cex <- 0.66
	}
	if (K==171){
		dot.colour = "pink"
		line.colour <- "red"
		dot.cex <- 0.5
	}
	if (K==353){
		dot.colour = "cornflowerblue"
		line.colour <- "blue"
		dot.cex <- 0.33
	}	
	
	for (base in 1:T){
			points(abs(1:T-base),value.matched.models[[K]][base,],col=dot.colour,pch=16,cex=dot.cex)
			x <- c(x,abs(1:T-base))
			y <- c(y,value.matched.models[[K]][base,])
	}
	
	mean <- NULL
	for (diff in 0:29){mean <- c(mean,mean(y[which(x==diff)]))}
	lines(0:29,mean,col=line.colour)
}


quartz(width=9,height=7)
plot(NA,NA,typ="l",xlab="Quarter lag length",ylim=c(0,1),xlim = c(1,T), ylab="Value share, matched characteristics")
dot.colours <- rep("black",353)
c("grey","pink","cornflowerblue")
line.colours <- c("black","red","blue")
for (K in c(36,171,353)){
	y <- NULL
	x <- NULL
	if (K==36){
		dot.colour = "grey"
		line.colour <- "black"
		dot.cex <- 0.66
	}
	if (K==171){
		dot.colour = "pink"
		line.colour <- "red"
		dot.cex <- 0.5
	}
	if (K==353){
		dot.colour = "cornflowerblue"
		line.colour <- "blue"
		dot.cex <- 0.33
	}	
	
	for (base in 1:T){
			points(abs(1:T-base),value.matched.characteristics[[K]][base,],col="grey",pch=16,cex=0.5)
			x <- c(x,abs(1:T-base))
			y <- c(y,value.matched.characteristics[[K]][base,])
	}
	
	mean <- NULL
	for (diff in 0:29){mean <- c(mean,mean(y[which(x==diff)]))}
	lines(0:29,mean,col="black")
}
dev.off()

graphics.off()








